Mean-field analysis of collapsing and exploding Bose-Einstein condensates 
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The dynamics of collapsing and exploding trapped Bose-Einstein condensates caused by a sudden 
switch of interactions from repulsive to attractive are studied by numerically integrating the Gross- 
Pitaevskii equation with atomic loss for an axially symmetric trap. We investigate the decay rate 
of condensates and the phenomena of bursts and jets of atoms, and compare our results with those 
of the experiments performed by E. A. Donley et al. [Nature 412, 295 (2001)]. Our study suggests 
that the condensate decay and the burst production is due to local intermittent implosions in the 
condensate, and that atomic clouds of bursts and jets are coherent. We also predict nonlinear 
pattern formation caused by the density instability of attractive condensates. 



I. INTRODUCTION 

Bose-Einstein condensates (BECs) of trapped atomic 
vapor have been realized in several atomic species [pi @, 0, 
[|, a, n|. One of the remarkable features of these systems 
is that both strength |7J] and sign ||, 0, [10| of interac- 
tions between atoms can be tuned by adjusting an exter- 
nal magnetic field near a Feshbach resonance [JJTj. This 
has opened up the possibility to study collapsing and ex- 
ploding BECs in a controllable manner, thereby offering 
a stringent test of the Gross-Pitaevskii mean-field the- 
orygj. 

A trapped BEC with attractive interactions may be 
formed [ 13[ when quantum pressure arising from Heisen- 
berg's uncertain principle counterbalances the attractive 
force between atoms. This is possible when the parame- 
ter defined by k = No\a\/do is below a critical value fc cr , 
where No is the number of BEC atoms, a the s-wave scat- 
tering length, and do the size of the ground-state wave 
function of a harmonic trap. When the parameter k ex- 
ceeds fc cr , attractive force dominates quantum pressure, 
causing BEC to collapse. In the experiments performed 
by a Rice group |1J, |l^, |l6[ , BEC atoms are continuously 
supplied from a supercooled thermal gas, and collapse 
and growth cycles of BEC have been observed |L6|. In 
the experiments performed at JILA fi 0, HO], in contrast, 
No is fixed and a(< 0) is decreased by using the Fesh- 
bach resonance so as to meet the condition k > fc cr . This 
technique enabled them to determine the value of k CI [|9| , 
and to observe exploding atomic ejection from collaps- 
ing BEC [g, |l0| , a phenomenon called "Bosenova" whose 
origin is currently under controversy JlT], [L8|, [jjj ^Cj . 

The Gross-Pitaevskii (GP) equation |0| has widely 
been used to study mean-field properties of BECs, offer- 
ing fairly good quantitative account of a rich variety of 
experiments for repulsive BECs plj . In the case of the 
attractive BECs, however, it is by no means clear to what 
extent mean-field theory is valid, for the attractive inter- 
action might enhance many-body quantum correlations. 
In fact, a measured value of k CI H is significa ntly smaller 
than that predicted by a mean-field theory (l3| |22|, the 
origin of the discrepancy being not yet understood. 

The recent quantitative measurements on collapsing 



and exploding BEC reported in Ref. |Hj have motivated 
us to investigate to what extent mean-field theory can 
explain the experimental observations. This is the main 
purpose of this paper. By numerically integrating the GP 
equation with atomic loss, we show that the phenomena 
reported in Ref. jig] , such as decay of condensates, ejec- 
tion and refocus of atomic bursts, and jet formation, are 
reproduced by our numerical simulations. In addition, 
we predict that various patterns in atomic density are 
formed in the course of collapse. 

This paper is organized as follows. Section II briefly 
reviews the GP equation with atomic loss due to inelas- 
tic collisions and describes the method of analysis of our 
simulations. Section III reports our results of numerical 



simulations using the GP equation, and compares them 
with the experimental data of Ref. [J10| . Section IV pro- 
vides the summary of this paper. 



II. THE GROSS-PITAEVSKII EQUATION WITH 
ATOMIC LOSS AND METHODS OF ANALYSIS 

The GP equation describes unitary time evolution of a 
macroscopic "wave function" -0, and conserves the total 
number of atoms J \ip\ 2 dr. In reality, however, atoms are 
lost from the trap due to the two-body dipolar and three- 
body recombination losses. These effects may be taken 
into account by incorporating in the GP equation the 
imaginary terms describing these inelastic processes [n7l : 
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where Vtrap is the trapping potential, and Ki and K3 
denote two-body dipolar and three-body recombination 
loss-rate coefficients, respectively. The imaginary terms 
in Eq. (|l|) are phenomenologically introduced in order to 
take account of the fact that the two-body and three- 
body losses are proportional to the square and cube of 
the atomic density: 
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We assume that the atoms and molecules produced by 
inelastic collisions escape from the trap without affecting 
the condensate. The constants K 2 and K$ include Bose 
statistical factors 1/2! and 1/3!, respectively, which are 
needed for BEC [E3[ , and describe the loss rate per atom. 

In the situations we consider the two-body loss can 
be ignored. The two-body and three-body loss rates are 
given by i?2 = K2\ip\ 2 an( i -^3 — ^3|'0| 4 : an d their ratio 
by i? 23 = R2/R3 = K 2 /(K 3 \i)\ 2 ). They have the relation 
R2R23 =K2/Ks, which is always < 0.1 s" 1 for 85 Rb ^ 
and 7 Li |2qj . When the density is low and R2 ~ R3, 
R2 ^ 0.1 s _1 . Then both two-body and three-body losses 
can be ignored since we shall consider the time scale of 
< 10 ms. When the density is high, the three-body loss 
dominates the two-body loss, and thus the two-body loss 
becomes unimportant. In our numerical simulations we 
shall therefore ignore the two-body loss term in Eq. (Q). 
All results presented below are not affected if the two- 
body loss is taken into account. 

The three-body loss, in contrast, plays a crucial role 
in determining the collapsing dynamics of BEC. When 
implosion occurs, the atomic density becomes extremely 
high, and so does the three-body recombination rate, 
until it stops the growth of the density. The maxi- 
mum density in the process of implosion is determined 
by K 3 and a 0, g ||. The values of K 3 far from the 
Feshbach resonance have been measured for 87 Rb p6[ , 
23 Na 27 J, and 7 Li [29], and they agree with theoretical 
predictions |2q, g9| 3?J within a factor of ten. Near the 
Feshbach resonance, however, complicated behaviors of 
K3 are predicted |H| , with no precise experimental data 
available J24|. 

We performed numerical integration of the GP equa- 
tion (n ) using a finite-difference method with the Crank- 
Nicholson scheme fl3[| . Since the peak density changes 
drastically during implosion |18, [l9|, we very carefully 
controlled the time step to avoid error propagation. Ini- 
tially we prepared the ground-state wave function by the 
method in Ref. |32J for an initial s-wave scattering length 
Oinit and an initial number of BEC atoms A/q. At t = 
the interaction is suddenly switched from ai n it > to 
^collapse < Oj inducing collapse of the condensate. We use 
the parameters of 85 Rb and assume the same trap geom- 
etry (i.e. radial frequency lu±/2tt = 17.5 Hz and axial 
frequency lo z /2tt = 6.8 Hz) as that used in Ref. pp| . 

In Ref. [[10| , atoms after the collapse are classified into 
three parts: remnant, burst, and missing atoms. The 
remnant BEC is a dense atomic cloud peaking around 
the center of the trap, and the burst is a dilute one that 
spreads broadly around the remnant BEC. In our simula- 
tions, we identify the remnant, burst, and missing atoms 
as follows. The one-dimensional density distributions de- 
fined by paxiaitz) = / \ip\ 2 dxdy and p Ia diai(r) = J \ip\ 2 dz 
show bimodal structures, one peaking at the center and 
the other spreading around it (an example is shown in 
Fig. 0b). We identify the former as the remnant and 
the latter as the burst, and determine the axial and ra- 
dial coordinates, Zb and r^, of the boundaries between 




FIG. 1: Time evolutions of the peak height of the wave 
function |?/>peak| (solid curve in arbitrary scale), the fraction 
Ntot/No of the total number of atoms remaining in the trap 
(dashed curve), and the fraction iV rcm nant /No of the remnant 
BEC atoms (circles) . The dotted curve is the best fit of the 
circles to Eq. (J4J) with i C oiia pS c = 3.6 ms and r dc cay = 3.7 
ms. At t — the s-wave scattering length is changed from 
«init = 7ao to a C oiiapsc = — 30ao , where ao is the Bohr radius. 
The initial number of BEC atoms is No = 15000, and the 
three-body recombination loss-rate coefficient is K3 = 2 x 
10 _28 cm 6 /s. 



the remnant BEC and the burst. We then calculate the 
number of atoms in the remnant BEC as 
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The burst atoms are defined by the ones outside the 
boundary. Because of ambiguities in defining the coordi- 
nates of the boundaries, N rcmnant and A^burst defined this 
way have uncertainties of about ±0.057Vo, while the total 
number of atoms in the trap iVtot = N rcmnant + A^urst is 
well defined. We define the number of missing atoms as 

^missing = Nq — N tot . 



III. RESULTS OF NUMERICAL SIMULATIONS 

A. Decay of the condensate and local intermittent 
implosions 

We first study the situation in which the scattering 
length is switched from aj n it = 7ao to a co iiapsc = — 30ao 
at t = 0, where ao is the Bohr radius. Figure l] shows 
time evolutions of the peak height iV'pcakl of the wave 
function, the total number of atoms in the trap AT to t, 
and the number of remnant BEC atoms iV remnan t, where 
we assume K 3 = 2x 10 _28 cm 6 /s. 

Our numerical simulations show that the condensate 
first contracts slowly with its peak height IV'peakl gradu- 
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FIG. 2: A snapshot of the column density of the implod- 
ing BEC taken at 2.5 ms after the s-wave scattering length is 
changed from ai n j t = to a co ii a psc = — 30ao. The initial num- 
ber of BEC atoms is iVo = 15000, and the loss-rate coefficient 
is K 3 = 2 x 10" 28 cm 6 /s. 



ally increasing. The total number of atoms remains con- 
stant during this process, since the recombination loss is 
negligible at such low densities. At t ~ 3.6 ms, implosion 
suddenly occurs in a very localized region (the size of 
the density spike is ~ 0.1 /mi while the size of the BEC 
cloud is several micrometers) for a very short period of 
time (~ 0.1 ms). Furthermore the implosion occurs not 
just once but many times intermittently for about 10 
ms fig , |l9fl . The three-body recombination loss promi- 
nently occurs during the implosion, since the atomic den- 
sity becomes extremely high. Several tens of atoms are 
lost in each intermittent implosion, resulting in a step- 
wise decrease of JV to t- The decay of JVremnant is due to 
both three-body recombination and atomic burst ejec- 
tion. We note that the behavior of N lemnan t, shown as 
the plots in Fig. [l], is very similar to the experimental 
result (Fig. lb of Ref. jHj). 

The implosions occur not only at the center of the 
trap but also at other locations on the trap axis, and 
more than one density spike is often seen simultaneously. 
A snapshot of the imploding process with ai n it = 0, 
^collapse = — 30ao, and No — 15000 is displayed in Fig. g, 
where the image is taken 2.5 ms after the switch of inter- 
actions. Two large spikes are seen on the trap axis. In 
an isotropic trap, on the other hand, implosions always 
occur one by one at the center of the trap. 
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FIG. 3: The collapse time t C oiiapsc at which implosions begin 
as a function of |a co ii apa e|/ao, where the initial number of BEC 
atoms is No = 6000. At t = the s-wave scattering length 
is switched from ainit = to a co iiapsc < 0. The dashed line 
indicates the critical ratio |a cr |/ao for iVo = 6000, below which 
the condensate does not collapse. 



Fig. EL icdiapsc = 3.6 ms, which agrees with the experi- 
mental result of 3.7 ms [M. Figure || shows the depen- 
dence of ^collapse on |a col ia P sc|/ao, where a in it = and 
Nq = 6000 [ p3[ . The plots are in good agreement with 
the experimental ones (Fig. 2 of Ref. [|l0| ) . For another 
set of parameters ainit = 89an and a co iiapsc = — 15arj, we 
obtain £ C oiiapsc = 15.9 ms, which is also consistent with 
the experimental finding mn . The collapse time increases 
for larger values of Oi n it , since the atomic cloud spreads 
more widely and therefore it takes longer time for the 
cloud to get to the center of the trap. 
Fitting 



Nremnant(*) = [No - iVremnant (oo)] e 

+A r remnant(oo) 



(4) 



to the plots in Fig. [0 (dotted curve) , we obtain the decay 
time constant Tdecay — 3.7 ms. For a- ln it — 0, and iVo = 
6000 and 15000, Td cca y is almost the same as the above 
one, and for a co n apso = — 60an it becomes Tdecay — 3.1 
ms which agrees reasonably well with the experimental 
finding of 2.8 ms p|. 



B. Collapse and decay times 

The total number of atoms of the system remains con- 
stant for some time after the switch of interactions, and 
suddenly it begins to decay. We call the time at which the 
sudden decay begins 'collapse time' £ C oiiapsc (see Fig. [l]). 
Since the collapse time is determined mainly by slow 
accumulation of atoms towards the center of the trap, 
^collapse only weakly depends on the value of K%. In 



C. Remnant, burst, and missing atoms 

Figure ^a shows the fractions of remnant, burst, and 
missing atoms after the implosion has finished for a; n it — 
and a C oiiapsc = — 30an- We note that AT reinrian t is much 
larger than N cr when Nq is large, as observed experi- 
mentally po| . To put it differently, once the condensate 
expands following the collapse, the critical density will 
not be reached even when iV remnan t > iV cr . This is be- 
cause the ratio luj_/lu z is irrational and therefore the axial 
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FIG. 4: The fractions of the remnant JV rGmnant /Ao (circles), 
burst NhuTBt/No (squares), and missing atoms A m i SS m g /Ao 
(triangles) after the collapse for (a) an axi-symmetric trap 
and (b) an isotropic trap, where No is the initial number of 
BEC atoms. The dotted curves show the fraction N cr /No of 
the critical number of atoms. The s-wave scattering length is 
switched from ainit = to a C oiia P sc = — 30ao with the loss-rate 
coefficient K 3 = 2 x KT 28 cm 6 /s. 



and radial refocuses of the burst atoms do not occur si- 
multaneously, with no further implosions occurring. 

This result presents a sharp contrast with that for an 
isotropic case shown in Fig. ||b, where the trap frequency 
is chosen to be the geometric mean (lu^luz) 1 ^ 3 . The num- 
bers of remnant atoms iV remnan t are always below A cr . 
This is because all collapsing atoms go to the center of 
the trap simultaneously, and therefore more implosions 
occur in an isotropic trap, which increases N m i ss [ ng and 
decreases iV remnant . Moreover, the implosions occur also 
when the burst atoms refocus, since they concentrate at 
the center of the trap. The data in Fig. Eh are taken be- 
fore the first refocus. Nevertheless, we note that N rcmnant 
is already below N„. 

We also note that in the axi-symmetric trap the frac- 
tions N rcmnant /N Q , N huist /N Q , and N misaing /N are al- 
most independent of No, particularly for TVq > 6000, 
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FIG. 5: (a) The integrated column density p(x,z) = 

J \ip\ 2 dy seen from the direction perpendicular to the trap axis 
and (b) the one-dimensional density distribution /5 ax i a i(z) = 
J \ip\ dxdy along the axial direction. The s-wave scattering 
length is switched from ajnit = to a co ii apsc = — 30ao with 
Ki = 2 x 10 _28 cm 6 /s, where the images are taken at t = 33.6 
ms. The images are smoothed in accordance with the exper- 
imental resolution (7 /im FWHM), and the central peaks are 
truncated. 



which is consistent with the experiments [|10fl . This is 
a consequence of the fact that the number of implosions 
occurring in the collapse is roughly proportional to iVo 
and that the numbers of the burst and missing atoms in 
each implosion are almost constant. 



D. Atomic bursts and "jets" 

A burst atom cloud is usually too broadly spread and 
hence too dilute to be seen. However, at every tt/loj_ 
(or tt/ll> z ) the cloud refocuses along the axial (or radial) 
direction and can be observed. Figure |5|a shows the col- 
umn density seen from the direction perpendicular to the 
trap axis and Fig. |5|b shows the one-dimensional distribu- 
tion yO ax iai(^) = / \iJA 2 dxdy, when the burst atoms focus 
along the trap axis (corresponding to Fig. 3 of Ref. Jul). 
The s-wave scattering length is switched from a; n it = 
to a co ii ap sc = — 30ao at t = and the image is taken at 
t = 33.6 ms ~ tcoiiapse + n/u>±. There are small peaks 
in the ridge of focus, which appear to correspond to the 
shoulders of the central peak in Fig. 3c of Ref. [pjj, Ac- 
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FIG. 6: The integrated column densities p(x,z) — J \ip\ 2 dy 
of condensates and jets seen from the radial direction. The 
wave function in Fig. H is expanded by a switch of the s-wave 
scattering length from a co iiapsc = — 30ao to a eX pand — lOOao- 
The size of each image is 48 x 90 /im. (a) The image is taken 
with i C xpand = 3.6 ms. (b) The sensitivity of imaging in (a) is 
increased, (c) The image is taken with i expan d = 5.5 ms and 
with the same sensitivity as (b). 



radial direction (the main jet at the center and two tiny 
jets on either side of it). The images (b) and (c), where 
the sensitivity of the imaging is increased, show that the 
jets are interference fringes. The parallel fringe pattern 
is characteristic of the interference between waves em- 
anating from two point sources |34| . In fact, there are 
two spikes in the atomic density that play the role of 
two point sources of matter waves in Fig. g. We find 
that the spacing between the fringes is proportional to 
^expand- This supports our interpretation, since the spac- 
ing between the fringes from two point sources is known 
to be given by ht cxpan( i/md J34| , where d is the distance 
between two point sources. In the experimental images in 
Ref. [jl(| (particularly in Figs. 5d-f), the jets seem to be 
ejected from the edge of the remnant condensate, where 
the density is too low to form the spikes, which also sup- 
ports our interpretation that the jets are the consequence 
of interference. Thus, the experimental observations of 
jets indicate that atoms expanding from spikes are coher- 
ent, suggesting that the burst atoms are also coherent. 



cording to our theory the origin of the atomic burst is 
a release of kinetic energy of local spikes in the atomic 
density 0,|l|,|l|. 

As regards the burst energy, it is at present difficult 
to compare our results with the experimental ones. One 
reason is that although the burst energy depends not 
only on a co n a p Se but also on K3 |l9| , experimental values 
of K$ as a function of a co n apsc are not available. The 
uncertainty in the boundaries between the remnant BEC 
and the burst and that in Gaussian fitting of the burst 
profile (see Fig. pp) give rise to large errors in numerically 
determining the burst "temperature" , which also make 
quantitative comparison between theory and experiment 
difficult. More work (both theoretical and experimental) 
needs to be done in order to clarify the situation. 

In the experiments [[To] , prolonged atomic clouds called 
"jets" are observed in the radial direction when the col- 
lapse is interrupted by switching the interaction from at- 
tractive to repulsive. The jets are distinguished from the 
bursts in that the energy of the former is much lower than 
the latter and that the direction of the jets is purely ra- 
dial. In Rcf. pi, the origin of the jets is considered to 
be the highly anisotropic spikes in the atomic density. 

We performed numerical simulations under situations 
similar to the experimental ones. The s-wave scattering 
length is switched from a lni t = to a co n ap sc = — 30ao 
at t = 0, and changed to a cxpan d = lOOao at t = 2.5 
ms (a snapshot before expansion is shown in Fig. g). 
Figure || shows the gray-scale images of the integrated 
column density p(x, z) — J\ip\ 2 dy after a Expand = 3.6 
ms expansion in (a) and (b), and i expa nd = 5.5 ms in 
(c). In the image (a), prolonged atomic clouds similar 
to the jets reported in Fig. 5 of Ref. Gof appear in the 



E. Pattern formation 



We found in Ref. [18J that various patterns in the 
atomic density are formed in the collapse processes 
caused by a sudden switch of interactions from repulsive 
to attractive, the origin of the pattern formation being 
attributed to the self- focusing effect of the attractive sys- 
tems. Here we predict that a similar pattern formation 
does occur in 85 Rb BEC for experimentally available pa- 
rameters. 

Figure [fl shows pattern formation in axi-symmetric 
traps, where iVo = 5 x 10 4 , and the s-wave scatter- 
ing length is switched from ai n ; t = 400ao to a co n apsc = 
— 310ao. The pancake-shape trap has the ratio Wz/wj. = 
Vo with the same geometric mean frequencies as the 
cigar-shape trap. In the cigar-shape trap the cylindri- 
cal shell structure is formed, and in the pancake-shape 
trap the layered structure is formed. 

Figure ^ shows time evolution of the integrated column 
density p(x,z) = J \ip\ 2 dy for an isotropic trap, where 
N a , a in it, and a co n apsc are the same as in Fig. 0. We see 
that the density fluctuations grow to form four concentric 
spherical shells at t — 7.7 ms, which move inwards and 
collapse one by one; then at t = 11 ms we see that a new 
shell is being formed 

The resolution of the imaging system in Ref. |1(| (7 
/an FWHM) is inadequate for observing the patterns in 
Figs. M and H (the spacing between the shells is ~ 1 /j,m). 
Expansion of BEC before imaging will blur out the pat- 
tern. In order to observe the pattern formation, there- 
fore, we need to improve the in situ imaging method, or 
to use larger do and \a\ to enlarge the pattern. 
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The validity of the mean-held GP equation is deter- 
mined by the gas parameter na 3 , and the depletion is 
given by ~ (na 3 ) 1 / 2 . When the implosion occurs, n|a| 3 
becomes ~ 10~ 3 at the peak density in our simulations, 
which indicates that the mean-field approximation is still 
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FIG. 7: Pattern formation in a cigar-shape trap (left panel) 
and a pancake-shape trap (right panel), where the column 
density is seen from the radial (upper images) and axial 
(lower images) directions. The initial number of atoms is 
No = 50000, and at t — the s-wave scattering length is 
switched from ainit = 400ao to a C oiia P se = — 310ao. The ratio 
is lu z /lu ± = 0.39 (V8) with u± = 17.5 Hz (9.03 Hz) in the 
cigar-shape (pancake-shape) trap, and the images are taken 
at t = 6.2 ms (t - 4.8 ms). 
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IV. SUMMARY 

We have studied the dynamics of collapsing and ex- 
ploding BECs by numerically solving the time-dependent 
GP equation with atomic loss (Of), and compared our re- 
sults with those of the experiments of Ref. |Hj . We find 
that mean-field theory with atomic loss can account for 
the following experimental findings: (i) It takes the sys- 
tem a certain time £ C oiiapse to undergo a sudden decrease 
in the number of atoms after the jump of a. (ii) The 
number of atoms in the condensate decays exponentially 
with a decay time constant of a few milliseconds, (iii) 
The burst atoms are ejected in the collapse process, and 
refocus after every half trap period, (iv) The fractions of 
remnant, burst, and missing atoms are almost indepen- 
dent of No, and the number of remnant atoms is much 
larger than the critical number iV cr for large Nq. (v) 
The jets are observed when the collapse is interrupted by 
jumping a to a positive value. 

We have found that these phenomena are attributed 
to a rapid sequence of local intermittent implosions, and 
provided a new interpretation of the jets, i.e., the highly 
anisotropy of the jets is due to the interference fringes. 
This suggests that the burst atom cloud is coherent. 



FIG. 8: Time evolution of the column density of BEC in 
an isotropic trap, where the trap frequency is lu/2tt — 12.8 
Hz. At t = the s-wave scattering length is switched from 
ainit = 400ao to a co n apsc = — 310ao. The initial number of 
atoms is No = 50000 and the loss-rate coefficient is Ks = 

26, 



8 x 10- 26 cm 6 /s 



valid at least qualitatively. 

Our results presented here suggest that the mean-field 
approximation can be used to describe the collapsing and 
exploding dynamics at least qualitatively. A more quan- 
titative comparison between experiments and numerical 
simulations might reveal effects beyond mean-field ap- 
proximation. This possibility merits further experimen- 
tal and theoretical study. 
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